% This file loads in existing figures and adjust infrastructure quality in
% rural areas

load conus
earthradius = almanac('earth','radius');

countries = {'Bolivia'};
for c=countries
    load([path_load_grids, cell2mat(c), '_grid_0.2_0.2_0.5_connected1_WP0_osm1.mat' ]);
    gridmap = country_graph.gridmap;
    graph_export = country_graph.graph_export;
    
    % First identify rural areas
    % http://documents1.worldbank.org/curated/en/417881468226751166/pdf/wps3634.pdf
    for j=1:length(gridmap) 
        gridmap(j).area = max(areaint( gridmap(j).Y, gridmap(j).X, earthradius)); % area is in sq-qm
        gridmap(j).popdens = gridmap(j).population/gridmap(j).area;
        gridmap(j).rural = gridmap(j).popdens<20;
    end
    
    % Share of population living in a rural area
    sum(cell2mat({gridmap(:).population}).*cell2mat({gridmap(:).rural}))/sum(cell2mat({gridmap(:).population}))
    
    % then edit infrastructure quality matrix 
    count = 0;
    for o=1:length(gridmap) 
        for d=1:length(gridmap) 
            if graph_export.avI(o, d)>1e-4
                if gridmap(o).rural==1 && gridmap(d).rural ==1
                    graph_export.avI(o, d) = 0.8*graph_export.avI(o, d);
                    count=count+1;
                elseif (gridmap(o).rural==0 && gridmap(d).rural ==1) || (gridmap(o).rural==1 && gridmap(d).rural ==0)
                    graph_export.avI(o, d) = 0.9*graph_export.avI(o, d);
                    count=count+1;
                end     
            end
        end
    end
    
    % add back to original graph
    country_graph.gridmap = gridmap;
    country_graph.graph_export = graph_export;
    
    % write a new file 
    save([path_load_grids, cell2mat(c), '_grid_0.2_0.2_0.5_connected1_WP0_osm1_rural.mat'], 'country_graph' );
end